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A survey is given on the applications of hydrodynamic model of nucleus-nucleus collisons, focusing 
especially on i) the resolution of hydrodynamic equations for arbitrary configurations, by using the 
smoothed-particle hydrodynamic approach; ii) effects of the event- by- event fluctuation of the initial 
conditions on the observables; iii) decoupling criteria; iv) analytical solutions; and others. 
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I. INTRODUCTION 



T^. Kffont of resonance width a,s a, fimction of tomnoratn^ drodynamic Model has been proposed by Landau [l| 

m 1953 as an improvement over the Fermi statistical 
model for the multiple particle production phenom- 
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ena in high-energy nuclear collisions. At that time, these 
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f crmi model ottered an ingenious insight into the mech- 
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anism of the high-energy nuclear collision processes and 
gave a prediction for the energy dependence of the mul- 
tiplicity, which was verified by the data, it was known 
that it had troubles in reproducing particle spectra and 
relative abundance of K over tt. This was because, in 
this model, particles were assumed to be emitted directly 
from the hot and dense, thermally equilibrated matter 
formed in high-energy nuclear collisions, which was sup- 
posed to be at rest, so that the model predicted isotropic 
momentum distribution which did not agree with the ob- 

Dorvpd [Spectra. Furthermore, because of high tempera- 
Transverse expansion on longitudinal scaling exT|an^ic^^ §^^ ^^^^^^^ p^^^^^^^ multiplicity ra- 
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tio dopondo d only on the isotopic-spin statistical weight. 

Shock formation and Ncumann-Rachtmeyer pseujjo v^|^o^jA^A|l 6^ 4/3 r^^-^ conclusion of the model was 
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also not in agreement with data. 

These problems were solved naturally by letting the 
hot and dense matter expand before particle emission 
takes place, reducing thus heavy-particle multiplicities, 
because of the Boltzmann factor, and giving at the same 
time alongated momentum spectra, due to a violent lon- 
gitudinal expansion caused by a large pressure gradient 
in the beam direction. A nice feature of this model is 
that, since the entropy is conserved in the ideal case Lan- 
dau studied, the energy dependence of the total particle 
multiplicity predicted by the Fermi model, and verified 
experimentally, is preserved. 

When accelerator data on multiparticle production be- 
gan to appear, first in pp collisions at CERN ISR, and 
later in pp collisions at SppS collider, Carruthers [Sj re- 
vived this Heretical Model in 1974, showing that sev- 
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eral aspects of those phenomena may be well understood 
within Hydrodynamic Model. When laboratory study of 
high-energy heavy-nucleus collisions started, Hydrody- 
namic Model became one of the essential tools for these 
investigations. 

According to Hydrodynamic Model, the description of 
high-energy nuclear collisions goes as follows. At the 
beginning, two Lorentz contracted (in the cm. frame) 
nuclei collide and it is assumed that, after a complex 
process involving microscopic collisions of nuclear con- 
stituents, a hot and dense matter is formed, which would 
be in local thermal equilibrium. The description of this 
initial thermalization process is out of the scope of hy- 
drodynamics. In hydrodynamics, we simply assume that 
the local thermal equilibrium is attained and these states 
of matter are specified by some appropriate initial con- 
ditions (IC) in terms of distributions of fluid velocity 
and thermodynamical quantities for a given time-like 
parameter. Then, it follows a hydrodynamical expan- 
sion, described by the conservation equations of energy- 
momentum, baryon number and other conserved num- 
bers, such as strangeness, isotopic spin, etc. 

d.T^" = , (1) 

d^iuBu'^) = , (2) 

d^insu^) = , (3) 

where 

T^^ = (e + p)u''u'' - pgf"" (4) 

is the energy-momentum tensor, n^, ng, e, p are, re- 
spectively, the baryon number density, the strangeness 
density, the energy density and the pressure, all of them 
given in the proper frame of reference of the fluid ele- 
ment, and is the four-velocity of the fluid. Moreover, 
we have to specify some equations of state (EoS), which 
depend on the nature of the hot matter produced. 

As the expansion proceedes, the fluid becomes cooler 
and cooler and more rarefied, occurring finally the de- 
coupling of the constituent particles, that is, they don't 
interact any more until their detection. However, long- 
lived resonances and other unstable particles may decay 
after this instant of time. The observable quantities such 
as dN/dy, da/dniT, < V2 >, ■ • ■ are then computed by 
using these decoupled or free particles. 

The main object of studies by using Hydrodynamic 
Model is to investigate, through comparison of its pre- 
dictions with data, properties of the matter formed dur- 
ing high-energy nuclear collisions, specified by the initial 
conditions, equations of state and freeze-out or decou- 
pling conditions. We emphasize that these properties 
are not known a priori. It should also be stressed that 
even the basic assumption of "local equilibrium" is not 
granted for a priori. We expect that experimental and 
theoretical studies of some appropriate observables may 



respond these questions. Therefore, it is fundamental to 
find what are these "most appropriate observables" . 

In this survey, we shall discuss some aspects of this 
model, by focusing mostly on those ones, like develop- 
ment of hydrodynamic code capable to treat problems 
with highly asymmetrical configurations, effects of the 
initial-condition fluctuaions and improvement of the de- 
scription of decoupling process. These are features which 
have been investigated and developped within the Sao 
Paulo - Rio de Janeiro Collaboration in the last ^ 15 
years. For a review of other aspects of recent develop- 
ments, see for instance Ref. IMSH- 

In the following, in the next Section, we discuss the 
initial conditions, by emphasizing the importance of the 
event-by-event fluctuations as shown by some event gen- 
erators. In Section IIIIL we describe several equations 
of state, usually employed in these studies. Section Hvl 
is devoted to the resolution of hydrodynamic equations. 
There, we begin describing some analytic solutions, turn 
to the variational formulation and, flnally, an applica- 
tion of this approach to develop a numerical code, using 
algorithm of smoothed-particle hydrodynamics. Then, 
we consider the decoupling mechanisms in Section [3 by 
stressing that, although the commonly used Cooper-Frye 
presciption 4] is convenient and can give many good re- 
sults, more realistic treatment of decoupling is needed in 
order to correctly extract information on the hot matter 
formed in the collision process. In Section IVII we give 
some results obtained with the methodology described 
here. Finally, conclusions are drawn in Section IVIII 



II. INITIAL CONDITIONS 

In usual hydrodynamic approach of high-energy nu- 
clear collisions, one customarily assumes some highly 
symmetric and smooth IC, parametrized in a conve- 
nient way, which would correspond to the mean distri- 
butions of hvdrodynamic variables averaged over several 
events 0, 01 0- However, our systems are not large, 
so large fluctuations varying from event to event are ex- 
pected, even under the same initial conditions of colliding 
objects, such as the incident energy and the impact pa- 
rameter of the nuclei. What are the effects of the event- 
by-event fluctuation of ICl Are they sizable? Do they 
depend on EoS? These are some questions which arise 
regarding this subject. 

As mentioned in the Introduction, IC are determined 
by a complex process involving microscopic collisions of 
nuclear constituents not accounted for by hydrodynamic 
model, so when we want to introduce fluctuations in the 
IC of a hydrodynamic system, we must go beyond the 
hydrodynamic degrees of freedom. Just to see whether 
such event-by-event fluctuations of IC give sizeable ef- 
fects, so merit a more detailed study, in Paiva et al. 
used the Interacting Gluon Model [13 (IGM) to generat- 
ing fluctuating IC and, using Khalatnikov 1-dimensional 
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solution showed that the rapidity distribution ob- 
tained by averaging over results starting from fluctuating 
IC is quite different from that obtained starting from the 
averaged IC. 

There are some other simulations, which try to incor- 
porate, in hydrodynamic computations, fluctuating IC 
given by more elaborate microscopic models: with a use 
of some event generator, e^. HIJING [l^, VNI [13, 
URASiMA [lil, NeX uS lla . or some effective theory 
such as string model [l^, perturbative QCD -I- satu- 
ration of produced partons [T^ or color glass conden- 
sate jlS,]. In principle, one could test each of these dif- 
ferent microscopic models, by connecting them to some 
hydrodynamic code and computing several observables to 
see which are the differences among them and which are 
more suitable for describing experimental data, provided 
the other ingredients of the hydrodynamic model are well 
known, that is not the case. Here, we shall instead discuss 
not the details of such models, but more or less model- 
independent consequences of such fluctuations. Anyhow, 
we have to adopt some microscopic model. In the fol- 
lowing, we shall mainly discuss the recent works of Sao 
Paulo-Rio de Janeiro Collaboration, using NcXuS event 
generator, coupled to hydrodynamic code SPheRIO^. 

NeXuS is a microscopic model based on the Regge- 
Gribov theory and can give, in the event-by-event basis, 
detailed space distributions of energy-momentum tensor, 
baryon-number, strangeness and charge densities, at a 
given initial time r = \jW—z^ 1 fm, for any given 
pair of incident nuclei or hadrons. One important point 
when we use a microscopic model to create a set of IC for 
hydrodynamics is that the energy-momentum tensor pro- 
duced by the microscopic model does not necessarily cor- 
respond to that of local equilibrium. For example, NeXuS 
generates, as its output, the energy-momentum tensor 
T'^'^{x) and the current densities of conserved quantum 
numbers, jg{x), jg{x) and jqix), where B, S and Q re- 
fer to baryon number, strangeness and electric charge, 
respectively. However, the four-velocities corresponding 
to these currents usually do not coincide, and more im- 
portantly, do not coincide with that of the frame where 
T'^" becomes diagonal. Furthermore, the space compo- 
nents of the diagonalized T^"^ are not necessarily identi- 
cal (anisotropic stress). These facts mean that the mat- 
ter is not in local equilibrium. In order to transform the 
energy-momentum tensor to that of the equilibrated one, 
we adopt the following procedure. First, following Lan- 
dau 19], we identify the normalized time- like eigenvector 
of T'^" as the four- velocity of the fluid and the eigenvalue 
as the energy density, 

T'tw" = eu". (5) 



^ Smoothed Particle hydrodynamic evolution of Relativistic 
heavy lOn collisions. 



Using this four-velocity, we calculate the proper baryon 
number density as 

ns = JB^tJ. (6) 

and analogously for the other densities. Once e and n's 
are obtained, all the other thermodynamical quantities 
are calculated using the equations of state. By this proce- 
dure we force the system into a local thermal equilibrium, 
conserving the proper energy density of the system. 

The IC thus generated are used as inputs for SPhe- 
RIO code. We show in Figs. Q and |21 an example of 
such a fluctuating event, produced by NeXuS event gen- 
erator, for central Au -I- Au collision at 130A GeV, 
compared with an average over 30 events. As can 
be seen, the energy-density distribution for a single 
event (left), at the mid-rapidity plane, presents sev- 
eral blobs of high-density matter, whereas in the aver- 
aged IC (right) the distribution is smoothed out, even 




X (fm) X (fm) 



FIG. 1: Examples of initial conditions for central Au-f Au 
collisions given by NeXus at mid-rapidity plane. The energy 
density is plotted in units of GeV/fm^. Left: one random 
event. Right: average over 30 random events (corresponding 
to the smooth initial conditions in the usual hydro approach) . 




FIG. 2: A different representation of the same IC shown in 
FIG. Q at mid-rapidity plane. The vertical axis represents 
the energy density in units of GeV/fm'^. 

though the number of events is only 30. Similar bumpy 
event structure was also shown in calculations with HI- 
JING |l2j |. So, the main question here is whether the av- 
erages over the observables computed starting from such 
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fluctuating IC, like the one at the left: panel ofFigs.^and 
121 are similar or sizeably different from the correspondent 
ones computed from averaged smooth IC like that at the 
right panel there, or symbolically, 



1 ^ 



N 



J = l 



i=i 

(< IC >^ f) , 



N 



(7) 



where < / > is the value of some relevant quantity /, 
obtained by averaging over N total number of events 
and (IC— > f)j is the value of the same quantity in 
the j-th event, with some event-dependent IC, whereas 
(< IC >^ f) represents the same quantity / given by 
the average IC. This is a crucial point in data analy- 
ses, because the left-hand side is closer to the data point 
experimentalists obtain, whereas the right-hand side is 
the quantity usually computed by theorists for that data 
point, in order to extract properties of the matter formed 
in nuclear collisions. If the bumpy structure shown by 
event generators effectively exists in experimental situa- 
tions, then how do these hot spots manifest themselves in 
the observables? 

A general conclusion one can draw about this question 
is that the total entropy of the system becomes always 
smaller when one takes such fluctuations into account, 
in comparison to the case without fluctuations, which 
means with average over the event-by-event fluctuating 
IC taken before the expansion. This can be seen by ob- 
serving that, in ideal hydrodynamics, both energy and 
entropy are conserved. Then, considering for simplicity 
an ideal gas so that Si = a{Ei)^/'^, with a =const.> 0, 
for each random event. 
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(9) 



is proportional to the entropy for each particle species, 
one would expect that also the multiplicity becomes, in 
general, smaller when one takes such fluctuations into 
account. 

Other possible manifestations of this inhomogeneity of 
IC that we can expect intuitively are: enhancement of 
high-pT components due to more violent expansion in the 
surface region, smaller HBT radii, due to concentrations 
of matter in small spots, azimuthal asymmetry even in 
central collisions. However, computations are needed to 
obtain quantitative conclusions, whether such discrepan- 
cies are meaningful or not. We shall discuss this question 
in Sec. ED 



III. EQUATIONS OF STATE 

As mentioned already, the basic assumption in hydro- 
dynamical models is the local thermal equilibrium. Once 
this condition is satisfled, all the thermodynamical rela- 
tions should be valid in each space-time point ^. Thus, 
the energy, pressure and temperature are given as func- 
tions of baryon number and entropy densities, specifying 
the properties of the matter. In this Section, we discuss 
how to obtain simple phenomelogical equations of state 
(EoS) for the hydrodynamical description of relativistic 
nuclear collisions 1231. 



A. Hadronic gas 

The strong interactions among hadrons are very com- 
plicated and difficult to be incorporated into the EoS for 
practical use. However, for very high energy, we may 
consider that the hadronic gas may be approximated as 
an ideal gas, although the degree of approximation can 
not be evaluated theoretially. The recent thermal model 
for the description of chemical abundances '2^1 show that 
such an approach can reproduce quite well the observed 
multiplicity ratios of produced hadrons. Here we assume 
that all the particles can be treated as quantum ideal 
gas, except for a correction due to the excluded volume. 
We also include a main part of observed resonances in 
Particle Data Tables. The inclusion of resonances can be 
considered as an effective way to consider the interactions 
among hadrons as explained later. 

First, we recall that, in a grand canonical ensemble for 
an ideal gas of quantum particles, the thermodynamical 
potential per volume (the pressure) is given by 



Here, < E > and < S' > in the left-hand sides mean the 
averaged energy and entropy over the fluctuating events, 
whereas the right-hand side of < > is the entropy cor- 
responding to the averaged initial conditions, with the 
averaged energy < E >. The linear terms of the expan- 
sion in AEi/ < E > are cancelled out when the summa- 
tion is performed. If one recalls that particle multiplicity 
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(27r)3 



(10) 



^ Recently it is suggested that the thermodynamical relations can 
be satisfied without having the thermal equilibrium in the sense 



of Boltzmann distribution I20tl21l . 
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were 9 — ±1 (+ for fermions, — for bosons), [3 = l/T is 
the inverse of the temperature T, fi the chemical poten- 
tial, g the degeneracy factor and e(fc) = ^/k?■ + m? with 
m the mass of the particle. The number density n and 
the energy density e can be obtained by the usual ther- 
modynamical relations, n = {dp/dfi)v,T , e = {dp/dl3)\ , 
where A = e^^ is the fugacity. The entropy density of 
the gas can be calculated as s = /3(p + e — ^n). 

For example, in Landau's model 1], the equation of 
state was taken as that of the massless pion gas. For 
bosons with m — fi = 0, Eq. (|1U|I can be integrated 
analytically to give 



p{T) ^ i^T\ 
' 90 



(11) 



and accordingly 



and 



IT 



(16) 



where K2 is a modified Bessel function. From these rela- 
tions, we can see immediately the usual ideal gas relation. 



(17) 



When the widths of the resonances are taken into ac- 
count, Eg. 1)13(1 must be modified. For an interacting gas, 
the power series expansion of the pressure in terms of 
fugacity. 



T}_^bniT)e 



13 fin 



(18) 
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g TT 
~30 



rj-i4: 



(12) 



For the pion gas, due to the isospin factor, we can take 
9 = 5. 

For more realistic equations of state, we should include 
all the resonance particles in the gas. Furthermore, we 
should also take into account more than one type of con- 
served quantum numbers, such as electric charge (equiva- 
lently the 3'"'' component of isospin) , baryon number and 
strangeness. In this case, the chemical potential must be 
written as /i = BfiB + Sfis + r^^^Vs where B, S, T^^^ are 
baryon, strangeness and the thrid component of isospin 
quantum numbers, respectively, and ^J.B,^J■s and ^3 are 
the corresponding chemical potentials. Thus, for a mix- 
ture of particles with these conserved quantum numbers, 
Eq. H1U|) should be generalized to 



(13) 



is known as the cluster expansion (which is intimately re- 
lated with the virial expansion) and 6„ are called "virial" 
coefficients. Here, p"^ is the pressure of the corresponding 
ideal gas and, roughly speaking, the index n in the sum 
represents the order of multiple particle interactions. For 
n = 2, the contribution to the pressure comes from the 2- 
body interactions. Beth and Uhlenbeck [2^ showed that 
the second virial coefficient can be expressed in terms of 
the scattering phase-shift of constituting particles. This 
approach was generalized to the relativistic Boltzmann 
gas by Dashen, Ma and Bernstein (2^] and the result for 
62 is 



dW W^K2{fiW) 



(19) 



where 5i{W) is the phase shift for the ^-th partial wave. 



where the sum refers to the particle species (including res- 

onances) and /i^ — BifiB + Sifis+T^ /is with Bi, Si, T> 
are the quantum numbers of the i-th particle species. We 
verify that the baryon number density of the mixture is 



ub 



dp 



E 



dpt (T, lij) 
dfiB 



(14) 



where n^'-* — [dpi{T , ^i) / d ^i) is the number density of 
the z-th particle species. 

Except for pions, most of hadrons and resonances can 
be well approximated by the Boltzmann limit. In this 
case, we have 



(15) 



Consider the case of a gas with mass M and suppose 
there exists a resonance in the two particle collision. 



M + M Mr . 



(20) 



When the resonance has a width F and spin J, only £ = J 
dominates the sum and 



2Mr-W' 
so that we have the Breit-Wigner formula, 

dw^^ ' 2 (A/fl-iy)2 + r2/4 



(21) 



(22) 



Therefore, the pressure of the system can be written as 



P = P"^ + PR. , 



where 



PR = 5-R 



47^3 



dW 



Wo 



W'^K2{f3W) 



, (23) 
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with 



9R = 25 + 1, 
For extremely narrow resonances, F — > , 



PR. 9R 



R)e' 



(24) 



which is exactly the pressure of the ideal relativistic 
Boltzmann gas made of resonances with mass . 
Eq. H23() suggests that, for more general case, the effect 
of resonance width can be obtained by a convolution of 
the normalized mass spectrum f{M) of the resonance, 
with the pressure p^'^{M) of the ideal gas of mass M, as 



PR 



dMRf{MR)f%M). 



B. 



Effect of resonance width as a function of 
temperature 
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FIG. 3: Correction factor F as function of temperature T for 
the resonance p. 



Equation H23|l shows that the effect of the resonance 
width on the pressure of the gas is temperature depen- 
dent. To see this effect, let us introduce the quantity, 



F{T,Mr) 



mIK2{Mr/t) 

Jwo (Mr 



W^K2{W/T) 



so that 



Pr=Pr='' y< FiT,MR,r). 



(25) 



(26) 



Another way to see the effect of the width, we may in- 
troduce the effective mass of the resonance Me/ / defined 



by 



MlffK^iPMeff) = MIK2{I3Mr)F{T,Mr). 



Using this effective mass, we can write the resonance 
pressure as 



PR = 9R- 



27r' 



|iiif2(/3Me//)e^'^«, (27) 



that is, as if the pressure of ideal particle with mass Mf,f f . 
In Figs. l3l4l51 andl^ we show the temperature dependence 
of F and Me/ / for two tipical cases, one for the resonance 
p (light and narrow width) and the resonance / (heavy 
and large width) resonances. As we see, the ideal gas ap- 
proximation deviates substancially for low temperature, 
especially for large width particles. The ideal gas approx- 
imation is only valid for T . 
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FIG. 4: Effective mass of the resonance p as function of tem- 
perature T. The dark area corresponds to the resonance 
width. 



Excluded-volume correction 



From the analysis of thermal modesl'23|, it became 
clear that the ideal gas description requires a modifica- 
tion to adjust the size of the system. The volume to 
fit the particle abundances is found to be too small. To 
avoid this problem, the correction due to the excluded 
volume effect, like a Van der Waals hard core correc- 
tion is introduced |26j |. According to this prescription, 
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FIG. 5: Correction factor F as function of temperature T* for 
the resonance /. 
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FIG. 6: Effective mass of the resonance / as function of 
temperature T. The dark area corresponds to the resonance 
width. 



Eq. is modified by the following coupled equations 



PHG{T,HB,fJ.S,f^3) = ^pf{T,|2^), (28) 

1=1 

lli = iii- ViPHG , (29) 



where as before fii = Bi^s + Sifis + ^3 is the chemical 
potential and Vi is the excluded volume of the i-th hadron 
species. The superscript id refers to the ideal gas case. 
The above equations constitute an implicit equation for 
Phg so that these two equations are solved iteratively to 
obtain phg for a given set of parameters, T,^B,^J'S and 



/i3. The number density of the i-th hadron is given by 



n. 



(30) 



where is the ideal gas expression of the particle den- 
sity for the z-th particle species. 



D. Gas of Quarks and Gluons 

The simplest way to introduce the phase of quarks and 
gluons in the equations of state is the use of the MIT Bag 
model. The effect of gluon and quark condensate in the 
physical vacuum is expressed as the enerdy density of the 
vacuum (or vacuum pressure). Thus the energy density 
and pressure of an ideal quark-gluon gas calculated in the 
QCD vacuum should be modified according to the rule, 



e 
P 



e - 
P 



B, 
B. 



where B is the vacuum pressure. Note that this vacuum 
pressure has an analogous property of the cosmological 
constant A of Einstein. Now, when we consider just the 
u and d quarks and neglect their masses, we have 



Pqgp 



9q 

67r2 



1 



-/i„ H 



60 



90 



(31) 



with 



9q 
9G 



2x2x3, 
2x8, 



the statisfical factors of quarks and gluons. For quarks, 
we have Hq = /LtB/3 . For /ljb = , we have 



TT 

Aud) _ 07 Y rpi n 



(32) 



or effectively gqgp = 37 . To include the strangeness and 
also charge conservation, we proceed in the same way as 
the hadronic gas and we have 

^ „ ,,-4t^4- 

Pqgp{T,^J.B,|J■S,^^3.)' 



91 

67r2 



1 4 2 77r'*T'* 
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+ pf(r,M,)+pf(T,~M.) 



where 



+ 

2x3 and 

Md = 



90 
1 



(33) 
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E. 



Construction of equations of state for the 
practical use 



The expressions Eqs. (|28I29|I or Eq. H33|) are, however, 
not convenient for the use in hydrodynamical calcula- 
tions. This is because the variables in such calculations 
are the conserved quantum numbers and the entropy den- 
sity and not the chemical potentials and the temperature. 
So, we need to invert the relations. 



to get 



ns = n3(T,^B,/i5,^3) , 

S = S (T,^S,^5,/Z3) , 



MB = fiB{nB,ns,n3, s) , 
iJ-s = \J^s{nB,ns,n-i,s) , 
M3 = fJ'3{nB,ns,n3,s) , 
T ^ T {nB,ns,n3,s) . 



(34) 
(35) 
(36) 
(37) 



(38) 
(39) 
(40) 
(41) 



However, this is a formidable task even numerically. We 
are thus forced to reduce the degrees of freedom for the 
practical application to hydrodynamics. For this pur- 
pose, we set the isospin and strangeness densities to null 
everywhere. That is, we impose the conditions, 



"-5 = 0, 
na = 0. 



(42) 
(43) 



These conditions together with Eqs. and deter- 
mine fis and fi3 as functions of T and fiB ■ Therefore, 
ub and s in Eqs. (|34|l and H37|l become now functions of 
two variables T and fiB , 



riB = nB{T,^B) , 
s = s{T,fiB) , 

which can be inverted numerically and give 

T = T{nB,s), 
fJ-B = HBinB,s). 



(44) 
(45) 



(46) 
(47) 



The above inversion process allows us to write any of 
thermodynamcial quantities as functions of ns and s, 
both in hadronic gas and quark-gluon plasma. Then, 
for a given pair of density parameters ns and s, we 
should determine which phase the physical system as- 
sumes. When two phases are in equlibrium, we must 
have m 



Phg{T,hb) ^pqgp{T,hb) 



(48) 



so that it determines the phase boundary line in (/is , T) 
plane and separates this plane into two domains. The 
domain where phg > Pqgp is the hadron gas phase and 
the other p^c < Pqgp is the QGP phase. These two do- 
mains are in contact on the phase boundary line Eq. (|48|l . 



However, the above two domains in (/xb,T) plane are 
mapped into two separated domains in the (ns, s) plane 
and there appeas a new third domain between them. 
That is, the phase boundary line (/iB,T) plane spreads 
into a domain in the (nB,s) plane. This domain is the 
mixed phase. In oder to determine thermodynamical 
quantities in this mixed phase as functions of (nB,s), 
we should introduce another criterion in addition to the 
phase boimdary condition of Gibbs. As mentioned above, 
any point (T, fii,) on the phase boundary line corresponds 
to two points in the (nB,s) plane, one for the hadron 
phase, {ng,s^) and other, (n^,s'3) for the QGP phase. 
In the mixed phase, any density of extensive quantity, 
should be a linear function of the qgp/hadron concentra- 
tion ratio. Thus for a given value of the baryon number 
density ub in the mixed phase, {ng,s^) and {n'^,s^) 
should satisfy 



{uB ~ nf ) + s". 



(49) 



From this equation, we determine the two points on the 
phase boundaries, {n^,s^) and {n^,s^). All the other 
extensive quantities, say e, can then be obtained as 



iff 



[ub ~ n'^) + e" . 



(50) 



Finally we can construct the equations of state in the 
whole region of (ns, s). The parameters of the final equa- 
tions of state are then, 

• Number of resonances included in the hadronic gas: 
Here, we take all the mesons with mass smaller than 
1.5 GeV, and baryons smaller than 2 GeV. Reso- 
nance widths are not included. 

• Quark masses: We may safely take m„ = rud — 0, 
but for the strange quark, we take = 120 MeV. 

• Size of excluded volume: In the example shown in 
the figures below, ~ (47rrQ/3) with ro = 0.7 fm 
for baryons and — Q for mesons. 

• Bag constant: We take B = 380 MeV/fm^. 

Fig.Qshows the line of constant temperature in (n^ , s) 
plane. Fig. |S1 shows the lines of constant entropy per 
baryon in (/is , T) plane. 

Note that the present equations of state have, by 
construction, the first order phase transition between 
hadronic gas and quark gluon plasma. Recent lattice cal- 
culations |28l | indicate that there exists a critical point so 
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Baryon Number Density (l/fm^) 

FIG. 7: Isotherms in (ns,s) plane. Values of the tempera- 
tures are indicated in the figure. Hadron gas, mixed phase 
and quark gluon plasma are clearly seen. 




% 400 800 1200 1600 2000 

FIG. 8: Constant entropy per nucleon curves in the {n, T) 
plane. The numbers indicate (s/ub) values. 



that the phase transition in the region starting from that 
point to zero /i b axis is not the first order but may be ei- 
ther a higher order transition, or a smooth cross over |[2^ . 
Effect of such possibihties should be investigated within 
the hydrodynamical models. 



IV. RESOLUTION OF HYDRODYNAMIC 
EQUATIONS 

In general, exact analytic resolution of relativistic hy- 
drodynamics is a difficult task due to the highly non- 
linear nature of these equations. So, usually one resorts 
to numerical computations. However, since the analyt- 
ical studies are more transparent, it would be useful to 
find analytical soluions, even though they correspond to 
highly ideal cases. Khalatnikov's one-dimensional an- 
alytical solution to Landau's initial conditions 
namely, ideal fluid at rest in a Lorentz contracted thin 
spatial region, gave rise to a new approach in high-energy 
physics. The boost-invariant solution 30] found 20 years 
later, is frequently utilized as the basis for estimations 
of initial energy densities in ultra-relativistic nucleus- 
nucleus collisions jsj . Below, we shall first describe fam- 
ilies of analytical solutions we obtained in collaboration 

with T. csorgo mm. 

Considering that it is not trivial to analytically solve 
the equations of hydrodynamics, it would be nice if we 
could develop a method to obtain an approximate but an- 
alytical solution of hydrodynamics. This has been done 
with variational formulation [sil] . although we had not 
develop it further by applying it to practical problems 
of high-energy nuclear collisions. However, we did ap- 
ply the variational method to develop a numerical code 
SPheRIO, based on the so-called smoothed-particle hy- 
drodynamics (SPH) tecnique, to high-energy nu- 
clear collisions which is flexible enough, capable to 
treat systems with configurations without any symmetry 
and also exploding in time. 

In the following, after presenting some analytical solu- 
tions, we shall describe in the subsection IIV Bl the vari- 
ational formulation of hydrodynamics, showing how it 
could used to get approximate solutions. Then, in the 
subsection llV CI we shall apply this method to adapt the 
SPH hydrodynamics for relativistic heavy-ion collisions. 



A. Analytic solutions 

After Landau's initial proposal the first analytic 
solution obtained is due to Khalatnikov i considering 
1-dimensional expansion of ideal gas. A more simpler 
boost-invariant solution has been obtained later 0^ and 
appHed to estimate the initial energy densities in ultra- 
relativistic nucleus- nucleus collisions (3l| . Both of these 
have been frequently used in the study of nuclear col- 
lisions, showing the usefulness of such simple analytical 
solutions. 

However, the boost-invariant solution has some short- 
comings: it is scale invariant, having a flat rapidity dis- 
tribution, corresponding to the extreme relativistic colli- 
sions, which has never been seen; ii) it contains no trans- 
verse flow. In we tried to overcome these short- 
comings. We started by assuming a simple equations 



10 



of state, corresponding to a gas containing massive con- 
served quanta, namely. 



e = mn 
p = nT, 



(51) 
(52) 



having two free parameters, m and k. Non-relativistic 
hydrodynamics of ideal gases corresponds to the limiting 
case of m >> T, << 1 and k = 3/2 . 

Then, we looked for similarity flows, i.e., 



So, this is a generalization of the one-dimensional scale- 
invariant solution, including a class of transverse flows. 

More recently [s^l, we extended these solutions still 
further, considering less symmetrical flows, but still keep- 
ing the same EoS H1I2|I and similarity flow, with constant 
velocity, as it appears in eas. H5t)l57|) . 



B. Variational formulation 



x(i)"' r(t)"' z{t)' 

F{s), 



f ( 



(53) 
(54) 



where X{t),Y{t), Z{t) are the scales of length, in three 
orthogonal directions, and f{x^) is any of the thermody- 
namical quantities, such as n{x^),T{x^),p{x^^), . . . , and 



X 

X2 



Y2 



z 

Z2 



(55) 



is a scaling variable. 

Using this parametrization for the one-dimensional ex- 
pansion, it can be easily verified, by direct substitution 
into eas. (|ll2|) with eas. (|4I51I52I53I and . that a family 
of solutions can be written as 



V 

n 
P 
T 



z/t = tanh?7, 
no(to/T)V(s), 
Po(^o/r)l+l/^ 

ro(io/r)i/« 



V(.)' 



(56) 



with po = n^T{) and where V(s) is an arbitrary non- 
negative function of s = z'^/{Z'^t^). The index stands 
for the initial values. Thus, this is a family of one- 
dimensional similarity flows, but it is not scale invariant 
and n and T are not constant for a constant t. 

Next, considering cylindrically symmetric flows, with 
boost invariance along z direction (collision axis), we 
could find the following family of solutions, for transverse 
flows: 

V = r/t , for |r| < t, 
n = no{TM/T)fV{s), 
P = Po(r.o/r)3+3/«, 



T = To(r.o/r))3A^ 



(57) 



Here, po = noTa and V(s) is an arbitrary non- negative 
function of s {— rl/lR^T^) in this case), with rt = 



As shown above, even for very simple equations of 
state, analytic solutions are limited to special cases. For 
realistic situations, even the equations of state themselves 
are available only in the form of numerical tables. There- 
fore, the numerical resources are essential for realistic 
studies of hydrodynamical behavior of ultra-relativistic 
coUisional processes. However, it is well-known that any 
numerical method for partial differencial equations re- 
quires highly sophisticated techniques to avoid numeri- 
cal instabilities, and usually it needs a very large scale 
computation, especially when we want to describe cor- 
rectly explosive processes such as relativistic heavy-ion 
collisions. However, as emphasized in the Introduction, 
in the hydrodynamic approach of high-energy nuclear col- 
lisions, its main ingredients, i.e., the equations of state, 
the initial conditions and the freezeout conditions are not 
quite well known. In such a situation, we actually don't 
need a very precise solution of the hydrodynamic equa- 
tions, but a general flow pattern which characterizes the 
final configuration of the system as a response to given 
equations of state, initial conditions and the decoupling 
procedure. So, we prefer a rather simple scheme of solv- 
ing the hydrodynamic equations, not unnecessarily too 
precise but robust enough to deal with any kind of ge- 
ometry. From this point of view, we stressed in Ref. |34| 
the advantage of a variational approach to relativistic 
hydrodynamics. 

Although not commonly found in general textbooks, 
the variational formulation of hydrodynamics has been 
studied by several authors j^, ^3 ■ -f^^f. [l^ , starting 
from the action 



/= / d^x{-e}, 



(58) 



where e is the proper energy density, the relativistic hy- 
drodynamics was derived from the variational principle. 

Here, we show its derivation, generalizing to include 
the rotational flow. To do this, we introduce the two vari- 
ables, the proper baryon density, n, and entropy density, 
s, which satisfy the conservation laws. 



d^inu^') = , 



(59) 



^x"^ ^ip- and 7?o = V -^o + ^'^^ '^^ = Vf^ - z'^ . where u'^ is the four velocity of the held, satisfying the 
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normalization, 



U^Uf^ = 1 



(60) 



The functional form of the energy density, 
e — e {n, s) 

specifies the thermodynamical properties of the fluid. 
The pressure, temperature and chemical potential are ob- 
tained by the usual thermodynamic relations, 



P = 
T = 
fj. = 



d[e{n, s)/n] 
dn 
de{n, s) 

ds 
de{n, s) 
dn 



and 



The hydrodynamical equations of motion for the fluid 
is given by the variational principle with respect to n, s 
and under constraints, Eq. H59|l and the normaliza- 
tion of the four- velocity, Ea. H6U|) . These constraints can 
be incorporated in the variational principle in terms of 
Lagrangian multipliers to write 



- -w{u^u^- l)] = 0, 



(61) 



where A, C and w are Lagrangian multipliers and arbi- 
trary functions of x. Equivalently, the fluid dynamics is 
given by the effective Lagrangian, 

C[^fY'''\n, s, u^, A, Cw) = - e{n, s) - nu'^d^j^X 

(62) 

where now all of n,s,u^,X,(^,w are independent varia- 
tional variables. 

The variations with respect to n, s and u'^ lead imme- 
diately to 



-^-u^9^A = 0, 



(63) 
(64) 
(65) 



Variations with respect to A, C and w give simply the 
constraints, eas. H59l) and (|60|l . Multiplying the both sides 
of Ea. l|B5|l by u^, and using eas. H60l63l64f) . we get 



w = nji + Ts 



(66) 



where p is the pressure. Ea. H66|) shows that w is the 
enthalpy density. Sustiuting back this w into Ea. H65|) 
and multiplying u^, we have 

wui_,u^ = -{nu^){d^X) - (su^)(9^C) • 



Taking the divergence and using the continuity equations 
d'' (nu^) = and 9'^(su^) = 0, we get 

d'iwu^u,) = -(nu.Wd^X) - {su,){d''d^O ■ (67) 
But 

{nu^){d''dfj_X) = nd^^u^d^X) ~ n{d''X){dfj,u^) 
= -ndf_,fi - n{d''X){df_,Ui,) 

and analogously 

{su^Wd^O = -sd^T - sid^'Oid^u,) , 
so that Ea. (|67l) becomes 

= di_,p + u''{df^Uy)w 

= d^p. (68) 

Here, we have used Ea. H65|) and the Gibbs-Duhcm rela- 
tion. 



and the property. 



dp — sdT + ndfj, , 



(69) 



Finally we arrive at the standard form of the relativistic 
hydrodynamic equation 



where 



= {p + e)u^u„ - g^^p 



(70) 



(71) 



is the usual energy-momentum tensor of the fluid. 

It is important to observe that, the effective La- 
grangian Eq. (|62|l evaluated in the proper conioving frame 
of the fluid motion is 



C 



(fluid) 
eff 



comoving 



'e{n,s) + fm + Ts = p, (72) 



which is nothing but the negative of thermodynamical 
potential of the system. 

Now, interesting application of this approach appears, 
if we can parametrize possible solutions of continuity 
equations H59|) in terms of a certain number of time- 
dependent parameters d(t) = {ai(t), i = 1, . . . , N}, such 
that 

/ , da(t) \ 
n = „|r,a(t),^j. 
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together with the velocity field, 

then the action, Eq. (|58|l . may be written as a time 
integral of an effective Lagrangian 

Leff(a{t),'^)=- j dve{n,s). (73) 

The constraint terms vanish for these ansatz. In this 
case, the equations of motion for the variables ai{t) are 
obtained as the Euler-Lagrange equations. This method 
could be applied to relativistic heavy-ion collisions, try- 
ing to describe them in a simple analytic and effective 
way in terms of few parameters. However, the general 
parametrization which solve exactly the continuity equa- 
tions is not easy. An approximate way to solve the conti- 
nuity equation is proposed in a numerical method called 
Smoothed Particle Hydrodynamics. 

C. Smoothed particle hydrodynamics 

The SPH algorithm was first introduced for astrophys- 
ical applications [3^ [s^ . In [s^ , we extended this nu- 
merical method to heavy-ion collisions by the use of the 
variational approach discussed in the preceding subsec- 
tion. 



1. SPH representation of densities 

The basic idea of the SPH method is to parametrize 
the continous density distribution of any extensive phys- 
ical quantity in terms of sum of base functions with finite 
support. This procedure introduces two types of approx- 
imations of different nature. To see this, let us suppose 
that A is the physical extensive quantity and a(r, t) the 
corresponding density distribution. We start with the 
identity, 

a(r,i) = j a(r',t)5(r-r')d^r'. 

Now, let us introduce the first approximation. Substitute 
the Dirac (5-function by a smooth, normalized function W 
with finite support, say, h, and transform the density a 
to a as 

a(r,t) ^S(r,i) = J a{r' ,t)W{r - r' ; h)d^r\ (74) 

where as mentioned, W is normalized, 

J W{r - r'; h)d^r' = 1 . 

and having the property of finite support, 

W{r - r'; h) ^ 0, for \r-r'\> h. 



At this stage, the new density a(r, t) describes the 
smoothed part of the original density a(r, i). From the 
Fourier transform, we can see that for this smoothed den- 
sity, the Fourier components with large wave numbers, 
corresponding to 



vanish rapidly. In other words, the kernel function W 
serves as the short wavelength cut-off filter. Physically, 
it is useful to introduce such a filter, since we very often 
want to eliminate very short scale part in order to extract 
the global feature of the dynamics of the system. In 60 's, 
similar idea has been used to smooth out the spectrum 
density of the nuclear shell model to extract the collective 
liquid drop potential by Strutinski |4fl| |. 

Now, we introduce the second approximation. This is 
rather to do with the practical aspect, that is, to reduce 
the degee of freedoms envolved in the calculation. We 
replace the integral Ea. (|74|l by a finite sum over finite 
discrete set of points, {rj, i = 1, .., N}. 

N 

a(r, t) aspH (r,t) = ^A,H^(r-r,;/i), (75) 

i 

where the weight Ai should be chosen appropriately to 
minimize the difference between a(r,t) and asp/f(r, i) 
everywhere. The above expression means that we are 
representing the continous density as sum of finite num- 
ber of unit distributions (kernel) carrying the quantity 
Ai. These unit density distributions are centered at the 
position r.; . 

Finally, the correspondence, 

N 

a(r, t) ^ asPH{r, t) ^'^(^ (76) 

i 

can be considered as an approximation ansatz for 
the density a(r, t) with finite number of parameters, 
{Ai, Ti,i = 1.., N} . Due to the normalization of the ker- 
nel W , we have 

N 

I asPH{r,t)d^r ^"^A,, 

•' i 

SO that we should choose 

N 

Y^A^ 

i 

as the total value of the quantity A of the system. 

2. Solution of continuity equations 

For the application in Hydrodynamics, we can use 
these parameters as the variational variables so that they 
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are depending on time. When we deal with one or more 
extensive quantities, we usually choose one conserved 
quantity as the reference density, say p, and represents 
it by the SPH form, choosing appropriately the weights 
{ui} to get 

N 

PSHP (r,i) = ^^,W^(r-r,(t);/i), (77) 



and take Vt constant in time. Other extensive quantity, 
say A, is calculated as in Ea. H76|l with weights. 



A,; = 



(78) 



The quantity {a/ p)i then represents the quantity A for 
the unit reference quantity p at the position r = ri{t) . 
Note that the time dependence of the density in Ea. H77|) 
comes from those of {ri(t)} . In this sense, Ea. H77() can 
be understood as if {Yi{ty\ are Lagrangian coordinates 
attached to small volumes of the order of /i^, with some 
conserved quantity, such as baryon number or entropy in 
the adiabatic expansion. From now on, we refer to these 
unit density distributions characterized by their coordi- 
nates {vi} as "SPH particles". From Eq.jTSJl, the quan- 
tity Ai can be interpreted as the part of A carried by the 
SPH particle i. 

The most powerful point of the above scheme of SPH 
representation is that we can solve the continuity equa- 
tion in a very simple manner. Suppose M is a conserved 
quantity. Then the corresponding density p should sat- 
isfy the continuity equation. 



dt 



V-(pv)=0, 



(79) 



where v is the velocity field. The SPH expression for the 
current j = pv is 



so that 



V • jspff (r, i) - 51 ^» l^^VWir - r,) . 



On the other hand, from Ea. lf77|l . 



dt 



dt 

dr,{t) 
' dt 



VV7(r-r,(0), 



By inspection, if we identify 

dr,{t) 

then we can see that Ea. (|79|) is automatically satisfied. 



For the application to the relativistic heavy ion colli- 
sions, we can take the entropy and baryon number as the 
basic conserved quantities. Then, their densities (in the 
space-fixed frame) are parametrized as 



N 



s*{r,t) = Y."^ W{r-r,{t)) 

i 

N 

n*{r,t) = W{r-r,{t)) 



(80) 
(81) 



where Vi and bg are the entropy and baryon number 
attached to the i-th "particle". The total entropy and 
baryon number are then given by 



S = j s*(r,t) = ^ j/i . 

■i 

N 

B ^ d^Y n*{r,t) = . 



(82) 
(83) 



The proper densities of entropy and baryon number are 
related with these space-fixed frame quantities as 



s 
n 



-1 * 
7 s , 

-1 * 
7 n , 



where 7 = u° is the Lorentz factor asscociated with the 
fluid velocity. For the hydrodynamical description of nu- 
clear and hadronic collisions at ultra-relativistic energies, 
we prefer to use the entropy than the baryon number as 
the reference conserved number to write the SPH rep- 
resentation of any other extensive quantities. This is 
because, the baryon number may become zero but the 
entropy density never vanishes in the physically interest- 
ing region within a hydrodynamical description. 



3. SPH action and SPH equations 

In the variational derivation of this method, the set of 
time-dependent variables {Yi,i — l,...,n} are taken as 
the variational degrees of freedom and their equations of 
motion are determined by minimizing the action for the 
hydrodynamic system. Thus, SPH may be considered as 
an effective description, in which the coordinates {vi(t)} 
associated with "particles" are the optimal dynamical pa- 
rameters which minimize the model action. We observe 
that or {hi\ are not dynamical variables and are 
determined by the initital conditions together with the 
constraints for the variational procedure. 

The effective Lagrangian, Eq. H73I) , is rewritten in SPH 
representation as 



LsPH{{r^Yi\) = - ^ Vi{els*), = - ^ 



(84) 
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where Ei is the "rest energy" of the i-th "particle" . Then, 
the equations of motion are obtained from the usual vari- 
ational procedure. This leads to the following coupled 
equations 



V,Ty(r,-rj;/i) = 0.(85) 



4- General coordinate system 

The variational procedure can readily be extended to 
coordinate systems with a non-Cartesian metric. The 
use of generalized coordinate systems is particularly im- 
portant when we consider realistic initial conditions for 
simulations of RHIC processes. As we know, in a rela- 
tivistic heavy-ion coUisional process, the initial state is 
a cold, quantum nuclear matter. Just after the colli- 
sion, the hadronic matter stays at a highly off-shell state 
and the materialization occurs only after ~ 1 fm/c in the 
proper time. Therefore, the local thermodynamical state 
would emerge for some local proper time and not for the 
global space-fixed time t. Thus, it is important to choose 
a convenient coordinate system for the description of the 
relativistic heavy-ion collisions. For example, one often 
uses the hyperbolic time and longitudinal coordinates to 
be described later. 

Let us consider a general coordinate system. 



(86) 



However, in order to unambiguously define the conserved 
quantity, we consider only the case that the time-like 
coordinate is orthogonal to the space-like coordinates. 



5m0 



0. 



(87) 



The action principle for the relativistic fluid motion can 
be written as 34] 



61 = -6 



-ge^O, 



together with the constraint for the conserved entropy 
current. 



1 



(89) 



iV^S-f) + -L ^ a, {V^SJV^) = , (90) 



where 



(91) 



and we use the notation. 



T = x*^ , "f — vP . 



The generalized gamma factor 7 is related to the velocity 
Va through UpW^ = 1, so that 



1 



500 - w ^ g w 



(92) 



where — g is the 3x3 space part of the metric tensor. 
That is 



.900 
-g 



(93) 



Let us now introduce the SPH representation. We may, 
for example, express the entropy density by the ansatz 



or by 



57 



^SPH 



J2^zW{r~n{T)), (94) 



as well. These two possibilities, besides others, are simply 
different ways to parametrize the variational ansatz in 
terms of a linear combination of given functions W{'r — 
fair)). The most important property of an ansatz should 
be that W satisfies the normalization condition imposed 
by the basic conserved quantity. Since the total entropy 
is expressed as 



i 

the normalization of W should be taken to be 

(ffW{f-f') = 1, 
for the parametrization Ea. l94|l and 

(fr^/^W{r-f') = 1, 



(96) 



(97) 



(98) 



for the parametrization Eq. H95() . In the usual SPH calcu- 
lations, it is not desirable to introduce in W the space- 
time dependence through its normalization condition. In 
this respect, the most natural way to introduce the SPH 
representation is Ea. l(M|l . With this choice, the SPH ac- 
tion is given by 



iSPH 



J dr J (fx E 

i 

j dT^u 



-9 SI 



W{r-n{T)) 



(99) 
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The variational principle leads to the following equation 
of motion, 



J2 "i^i 

3 

li Si 
Vi 




where 



li Vl 



(Vgoo 

p + e 



(100) 



(101) 



and the operator V is just the simple derivative operator 
with respect to the coordinate variable in use. 

For ultrarelativistic heavy-ion collisions, a useful set of 
variables is 



— tanh — 
2 t- z 



(102) 
(103) 

(104) 



As mentioned above, the initial conditions for RHIC pro- 
cesses are specified in terms of the proper time rather 
than of the coordinate time t. The variable r is not ex- 
actly the physical proper time of the matter, but for the 
initial times it may approximate the proper time. 

The metric tensor for this coordinate system is given 

by 



500 



Since the metric is space independent, we can use the 
parametrization 




where 

and W is normalized as 



/"OO 

47r / q^dq W{q) = 1 , 
"'0 



The SPH equation becomes 



ar T 



where the r\ component of the momentum is related to 
the velocity d-q/dr as 



p + £ \ drj 



* s J dr 
whereas in the transverse direction, we have 

f p + e\ drx 

The Lorentz factor is given by 

1 

7 



1 — — T^v'^j 



Landau Model 



In order to show the efficiency of the method and also 
to show the corrrect choice of the coordinate system, let 
us investigate the Landau model in the SPH scheme, us- 
ing the ordinaly Cartesian coodinates and 77 — r cood- 
inates. Since the analytical solution is known, we can 
compare the numerical solutions to it. 

We thus solve the hydrodynamical evolution of a sys- 
tem of one-dimensional relativistic massless baryon-free 
gas initially at rest. The equation of state of a relativistic 
massless boson gas is 



where 



C 



( 



V 1287r2 



1/3 



To apply the SPH method, we introduce the discrete one- 
dimensional space variable Xi{t),i = 1, .., n (and similarly 
for r] — T coodinates) . The relation between the momen- 
tum and velocity is then 



(105) 



where, in this case, v can be solved analytically with 
respect to tt. 

In Figs.Ela and b, we show the results of our SPH cal- 
culation together with the exact solution 0. In these 
examples, we took only 100 particles with equally spaced 
Xi (or r]i). As we see from this example, in spite of rather 
small number of particles, the SPH solution is quite sat- 
isfactory. In particular, when we use the 77 — r coordi- 
nates with an appropriate distribution of i^'^s (Fig. |5l5), 
an excellent agreement with the analytical solution can 
be obtained. The computation time needed to get these 
solutions is even less than that needed to numerically 
evaluate the analytical solution. 



16 



0.4 




8 10 12 



t; 




FIG. 9: a) (above) Entropy profiles of tfie Landau model in 
Cartesian coordinate for different times. The exact results 
are given by the broken curves. The SPH solution is shown 
by the full curves, b) (below) Temperature profiles of the 
Landau model in the hyperbolic coordinate system (see text), 
for different time r. The SPH calculation is represented by 
the circles, and the exact result by the broken curves. 



6. Transverse expansion on longitudinal scaling expansion 

As a further test, closer to a realistic situation than 
that of Figs. El we calculated the transverse expansion 
of a cylindrically symmetric homogeneous massless pion 
gas, undergoing a longitudinal scaling expansion, and ini- 
tially at rest in transverse directions. Such a problem 
has been discussed by several authors as a useful base to 
understand the transverse expansion. In Fig. 1101 we com- 
pare our results (a full 31? calculation without assuming 
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FIG. 10; Temperature profile of a cilyndrically symmetric flow 
with longitudinally scaling expansion, shown as a function 
of r = ^ -^-y^ . The SPH results at = (circles) are 
compared with the numerical solution obtained by a space- 
fixed grid method. The SPH calculation has been perfomed 
in fuU 3D. 



cylindrical symmetry) with (2-1-1) numerical results, ob- 
tained by the use of the method of characteristics [4ll |. 
In this example, we used also 50 x 50 x 50 particles. The 
result is quite satisfactory. If we decrease the accuracy 
by 10%, we can reduce the particle number almost by 
one order of magnitude. 



7. Shock formation and N eumann-Richtmeyer pseudo 
viscosity 

As seen in the previous examples, our entropy-based 
relativistic SPH method works quite well for the adia- 
batic dynamics of the massless pion gas. However, for the 
application to realistic problems, it is fundamental to see 
how this scheme works for non-adiabatic cases, too. This 
is because, whenever a piece of fluid matter flows into an- 
other region of the fluid with a speed exceeding the sound 
velocity of the fluid, there appears a shock wave, and this 
is essencially a nonadiabatic process. Thus, except for a 
really quasi-static dynamics, there should be an entropy 
production mechanism. This becomes especially impor- 
tant in a domain close to the phase transition region, 
because there the velocity of sound tends to zero. In the 
following, we study some examples of one dimensional 
shock problems in the scheme of the SPH methods. 

The shock front manifests as a discontinuity in thermo- 
dynamical quantities in a hydrodynamic solution. Math- 
ematically speaking, the shock front should be treated 
as a boundary connecting two distinct hydrodynamic so- 
lutions. The smoothed particle ansatz excludes such a 
possibility from the beginning. Since short-wavelength 
excitation modes do not exist in the SPH ansatz. the 
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energy and momentum conservation required by the hy- 
drodynamics results in very rapidly oscillating motion of 
each SPH particle. Such a situation occurs, for example, 
when a very high-energy density gas is released into a 
low density region. This kind of shock, for the case of 
a baryon gas, is discussed in and also, in the SPH 
context, in |4J|. Here, we appliy our entropy-based SPH 
approach to the massless pion gas. 

Fig. 111! gives the typical behavior of SPH solution for 
such a situation, if entropy production is not taken into 
account. As discussed above, there appear in fact rapid 
oscillations in thermodynamical quantities just behind 
the shock front. Actually, such oscillations always appear 
in any numerical approach if entropy production is not 
included. 
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FIG. 11: Shock wave formation in one-dimensional pion gas, 
calculated with SPH. No viscosity is used. 

In order to avoid these unphysical oscillations, von 
Neuman and Richtmeyer 42] introduced the concept of 
pseudoviscosity. The idea is to set the dissipative pres- 
sure where the shock wave discontinuity is present. To 
do this, Neuman and Richtmeyer proposed to replace the 
pressure by 

where Q is the pseudoviscosity and they took the follow- 
ing ansatz. 



{aAx)^p {p/p)\ 
0, 



p > 
p < 



the space grid size and a is a constant of the order of 
unity. In order to generalize the above pseudoviscosity 
for relativistic SPH case, we replace the quantity p/ p 
by —9 = —df^u^ and Ax by h, where h is as before the 
width of the smoothing kernel W. More precisely, we 
take the following form which is a slightly modified ex- 
pression suggested by Ref. [i^. 



Q 



p 



-ah0 
0, 



<0, 

e>o. 



(106) 



where 



V dt 



(107) 



Actually, Q is equivalent to the bulk viscosity and 
therefore there is no heat flow associated with it. What 
this artificial viscosity does is to convert the collective 
flow energy into the microscopic thermal energy. As a 
consequence, the total energy, that is, the sum of the 
collective flow energy and the internal thermal energy 
is still conserved. In order to incorporate the internal 
energy conservation in the SPH scheme, we substitute 
all the pressure pi by pi + Qi, and we add the following 
equation for the entropy production. 



1 dui 
Vj dt 



Ts* 



(108) 



Fig-Elis the solution of the same problem as in Fig. 1111 

but with the entropy production taken into account. In 
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The above formula is for nonrelativistic one-dimensional 
hydrodynamics. Here, p is the mass density, Ax is 



FIG. 12: After the introduction of the Q term in the SPH 
calculation. 
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this calculation, the parameters have been chosen as 



a 



2, /3 



and h — Q.5fm for 1000 SPH particles. As we see, the 
rapid oscillations have been smoothed out (and in turn, 
the numerical calculation became much more efficient). 

It is known that the energy- and momentum-flux con- 
servations through a sock front relate the ratio S2/S1 of 
entropy densities after and before the shock to the veloc- 
ity Vs of the shock front as (Hugoniot-Rankine relation) 



Si 33/4'''' (l_«2)5/4 



(109) 



In Fig. ^1 we show the velocity of the shock front ob- 
tained in our SPH calculations as function of the en- 
tropy ratio(dots). Each point corresponds to the dif- 
ferent initial condition. They are compared with the 
Hugoniot-Rankine relation Eq. (|109|l (curve). The accor- 
dance shows that our SPH calculation reproduces faith- 
fully the conservation of kinetic energy and momentum 
of the flow through the shock front. 




1 2 3 4 5 6 7 8 
Ratio of entropy densities 

FIG. 13: Test of Hugoniot-Rankine relation. The circles are 
results of SPH calculations with different initial conditions, 
and the full curve is Eq. 111091 . 

In the usual hydrodynamic computations using space 
grids, the symmetry of the problem is often a crucial 
factor to perform a calculation of reasonable size. The 
SPH method cures this aspect and furnishes a robust 
algorithm particularly appropriate to the description of 
processes where rapid expansions of the fluid should be 
treated. The equations of motion are derived by a varia- 
tional procedure from the SPH model action with respect 
to the Lagrangian comoving coordinates. This guaran- 
tees that the method furnishes the maximal efficiency for 



a given number of degrees of freedom, keeping strictly the 
energy and momentum conservation. For this reason, so- 
lutions can be obtained with a very reasonable precision, 
with a relatively small number of SPH particles. This 
is the basic advantage of the present method, when we 
analyze the event-by-event dynamics of the relativistic 
heavy-ion collisions. 

On the other hand, the precision of this method in- 
creases rather slowly with the number of SPH particles. 
Therefore, a relatively large number of particles is re- 
quired if one wants a very precise numerical solution. 
However, for the application to the RHIC physics, we 
may need rather crude precision especially if we consider 
the dubious validity of the rigorous hydrodynamics. For a 
calculation with typically 10% errors, the SPH algorithm 
presented here furnishes a very efficient tool to study the 
flow phenomena in the RHIC physics. 

A fundamental difficulty of the relativistic hydrody- 
namics for viscous fluid l46l | is that the dissipation 
term causes an intrinsic instability to the system. This 
instability basically comes from the fact that the dissi- 
pation term contains 6 = d^u^ (see Eqs. (I106ll08|l '). so 
that it necessarily introduces the third time-derivative 
into the equation. This means that we have to specify, 
at least, a part of the acceleration as the initial condi- 
tion. Even we specify the initial acceleration, the require- 
ment of the internal self-consistency among the equations 
above leads to intrinsically unstable solutions. Israel pro- 
posed 0, 01 to cure these difficulties by introducing 
higher-order thermodynamics with respect to deviations 
from the equilibrium. Recently this "second order ther- 
modynamics" formalism was discussed in the context of 
Bjorken type solution ^^l- the examples presented in 
the present paper, we did not address this question and 
simply estimated the quantity 9 from the quantities one 
time step before. In practice, this will cause no numer- 
ical instability and the behavior of the solution is quite 
satisfactory. 

In spite of the above conceptual difficulties when nona- 
diabatic process is involved, the SPH approach has a nice 
feature as its flexibility, allowing the treatment of prob- 
lems with initial conditions without any symmetry, as 
happens in small systems as those resulting in relativistic 
nuclear collisions. It should be stressed that, due to the 
use of Lagrangian coordinates, the method is most suit- 
able for explosive processes like the relativistic heavy ion 
collisions. (Lagrangian coordinates have been used for 
treatment of relativistic nuclear collisions also by Non- 
aka, Honda and Muroya 0|). Furthermore, the varia- 
tioal approach guarantees that the SPH equations H85|) 
give the optimal description of motions for a given total 
number of "particles" {ri(t)}, which are our parameters. 
In this approach, no numerical instabilities will occur, 
since the whole system is a Lagrangian system. A nu- 
merical code, called SPheRIO has been developped by 
us on the basis of this algorithm. We shall discuss, in 
Sec lVIl some results obtained using this code. 
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V. DECOUPLING CRITERIA 

A. Cooper-Frye prescription 

As mentioned in the Introduction, the decouphng pro- 
cess is customarily described using the Cooper-Frye pre- 
scription PI , which gives the invariant momentum distri- 
bution as 



E- 



^(PN 



I daf.pf'fix^p). 



(110) 



This description of decouphng introduces a sharp freeze- 
out hypersurface a , usually characterized by a constant 
temperature Ty.o. . Before crossing it, particles have a 
hydrodynamical behavior and, when they cross it sud- 
denly decouple, free-streaming toward the detectors, 
keeping memory of the conditions (flow, temperature) 
of where and when they crossed the three dimensional 
surface. 

In SPH representation, we write 



E 



dp^ 



E 



(111) 



where the summation is over all the SPH particles, which 
should be taken where they cross the hyper-surface T = 
Tf o. and rij^ is the normal to this hyper-surface. 

Another often used procedure is to take such a freeze- 
out temperature not only constant for a given energy but 
also energy- independent. 

Though operationally simple, and actually useful for 
obtaining a nice comprehension of several aspects of the 
phenomena, such a concept of sharp freezeout hypersur- 
face and also of a constant freezeout temperature are 
clearly highly idealized when applied to finite-volume 
and finite-lifetime systems as those formed in high-energy 
heavy-ion collisions. 



Here, the energy dependence appears as a consequence 
of the increase in the initial energy density, which im- 
plies longer expansion time, so larger size L of the system 
(both longitudinally and tranversally) at the moment of 
decoupling, requiring lower density, so lower decoupling 
temperature, too. Figure H"!] shows the comparison made 
in [48j of an estimate of the incident-energy dependence 
of Ty.o. J using Landau's criterion mentioned above, with 
Tf o. obtained in a data analysis of tt and K transverse- 
momentum spectra in pp (and pp) collisions, in terms of 
a hydrodynamic parametrization of transverse velocity 
distribution and temperature. The data were taken from 
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FIG. 14; Energy dependence of T/.o. in pp (and pp) collisions. 
The solid line is the estimate as explained in the text. The 
points were obtained from a data analysis |43 |. 

An independent estimate of Ty.o.(s) for pp was made 
by Navarra et al. js^l, based on somewhat different but 
related freeze-out criterion 



B. Finite-size effect 

Before going further, let us for a moment assume that 
such a freezeout temperature is meaningful. At least, as 
an average temperature, it should exist. Then, how can 
we estimate it based on the properties of the system? 
A simple and natural criterion has already been given 
by Landau [J, by which a particle decouples when its 
mean free-path £ in the medium becomes larger than the 
system size L, 

i>L. (112) 

This means that 1/ o. is not an intrinsic thermodynamic 
property of the fluid, but it depends also on the size of 
the system. In fi^, we applied this idea to estimate Tj o. 
as function of the incident energy both for pp (pp) and 
nucleus-nucleus collisions, obtaining approximately 

Tf.o.-{V~^)-'/''- (113) 



'^hydro ^ '^col 5 

where Thydro and TcoI are, respectively, characteristic time 
for hydrodynamic expansion and particle collision, ob- 
taining a similar result. 

When Au -f An collisions data began to appear, Xu 
and Kaneta reported that ^ seems to decrease 
with yfs also in nuclear collisions at high energies, and 
we could verify that the reported results are consistent 
with TjFq. ~ (-^s)"^/^^. It is interesting that, more re- 
cently I54| . in trying to understand the experimentally 
observed anomalous behavior of the inverse slope param- 
eter T* of kaon transverse-momentum spectra in central 
Pb -f Pb (Au -f Au) collisions ^5^], we could succeed 
to obtain the increase in T* when going from SPS to 
RHIC domain only with decreasing freezeout tempera- 
ture Tf,o.{s) with increasing as shown in Fig. 1151 and 
Table HJ Here, SPheRIO code has been used with aver- 
aged IC. This is because, as will be shown in Sec. I VII 
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FIG. 15: Energy dependence of the inverse slope parameter 
for K'^ in central Pb + Pb (Au + Au) collisions. Here, av- 
eraged IC were used. Stars correspond to r/o=155 MeV (see 
Table EJ. Data are from js^ 

TABLE I: Freezeout temperature Tfo , averaged transverse 
velocity vt in —0.5 < y < 0.5 and the inverse slope param- 
eter T* for K'^ production in central Pb -I- Pb (Au -I- Au) 
collisions, obtained using SPheRlO with averaged NeXuS IC. 
To and eo are the initial values at the midpoint of the fluid. 
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Figs. El 12011^ Ell the slope parameter is not sensitive 
to IC fluctuations. 



C. Continuous emission 

The decoupUng temperature discussed in the preceding 
paragraphs should actually be interpreted as an average 
value of such a temperature. When applied to finite- 



volume and finite-lifetime systems as those formed in 
high-energy heavy-ion collisions, a sharply defined freeze- 
out hypersurface is clearly too idealized and for a more 
precise analysis of data, one needs a more elaborate de- 
scription of the decoupling process. In we intro- 
duced a method that we call Continuous Emission (CE) 
and, as compared to the usual Cooper-Frye one, we be- 
lieve closer to what happens in the actual collisions. 
The essential point is the introduction of momentum- 
dependent escape probability 



'P{x,p) — exp 



pavdr' 



(114) 



of a particle from a space-time point x without collision 
in the medium, so that the emission may occur from any 
point of the fiuid and at any time, according to this prob- 
ability. This means that we are interpreting probabilisti- 
cally the Landau condition, (|112|) . as should be and also 
giving the system size L a more precise meaning, namely, 
the quantity of matter the escaping particle encounters 
in his trajectory. The integral above is evaluated in the 
proper frame of the particle. Then, the distribution func- 
tion /(x, k) of the expanding system has two components, 
one representing the portion of the fiuid already free and 
another corresponding to the part still interacting, i.e., 

f{x,p) = ffree{x,p) + ftnt{x,p) . (115) 

We may write the free portion as 

ffreeix,p)^Vfix,p) . (116) 

The inclusive one-particle distribution, for instance, is 
then written as 



E- 



dp-^ 



daf,p'^ffree{xo,p) 
d^xdf,[p^ffree{x,p)], 



(117) 



where the surface term corresponds to particles already 
free at the initial time. 

As particles can be emitted in different stages of the 
fluid expansion, it is natural that the appearance of the 
observables in CE becomes different from that of the 
usual Cooper-Frye prescription at a constant tempera- 
ture 3. In general, in CE, the large-pr particles are 
mainly emitted at early times when the fiuid is hot and 
mostly from its surface, whereas the small-pT compo- 
nents are emitted later when the fluid is cooler and from 
a larger spatial domain. We shall discuss, in Sec. I VII the 
manifestation of this effect in two-pion interferometry. 

More detailed description of this method and several of 
its predictions are accounted for by Grassi in a separate 
paper of this issue . 



VI. APPLICATIONS 

In the following, we shall present some applications 
of NeXuS-fSPheRIO code El US HI, which has been 



21 



constructed by coupling the NeXuS event generator to 
SPheRIO code, described in Subsection IIVCI 

As mentioned in the Introduction, the main ingredi- 
ents of hydrodynamic calculations are IC, EoS and the 
decoupling procedure. None of these are well known at 
the present moment. 

Since we are mainly concerned with the effects of IC 
fluctuations and consequences of different decoupling de- 
scriptions, here we just take the IC created by NeXuS 
event generator and use the equations of state described 
in Section IIII El The strangeness conservation has not 
been taken into account. As for the decoupling proce- 
dure, we adopted both the conventional sharp freeze-out 
prescription and the continuous emission description, be- 
cause one of our purpose is to see the differences resulting 
from these two descriptions. 

We do not expect that these options will reproduce all 
the experimental data. In fact, we found that the use of 
the same version of NeXuS code for high- and low-energy 
nuclear collisions caused some discrepancies in reproduc- 
ing the rapidity distribution of charged particle. There- 
fore, we have introduced additional parameters to ajust 
at least the overall rapidity distributions. 



A. Effects of fluctuating initial conditions 

In Sec. inj we stressed that, due to the finite size of 

Rapidity Distributions (Pb+Pb, 17. 3A GeV) 
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FIG. 16: Rapidity distributions as indicated, with Tfo ~ 
140 MeV. The solid lines represent the averages over the fluc- 
tuating distributions (with dispersions), whereas the dashed 
lines are results with the averaged initial conditions. The data 
are from NA49j62i|. 

the systems, large fluctuations are expected in the initial 
stage of actual nuclear collisions, even for a fixed impact 
parameter, and that IC generated by realistic event gen- 
erators do show such effects 0, ^| . Let us see in this 
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FIG. 17: Centrality dependence of rapidity distributions, for 
Pb + Pb collisions at ^5 = 17.3^1 GeV (SPS). Here, "with 
fluctuation" (the solid lines) means results with NeXuS fluctu- 
ating IC and "without fluctuation" (dashed lines) those with 
the averaged IC. 



Section what are the effects of fluctuating IC on some of 
the observables. 

In T^l, by solving the hydrodynamic equations with 
longitudinal boost-invariance, it has been shown that the 
bumpy IC i) develop azimuthally asymmetric flows, even 
for central collisions, because there is no symmetry in 
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FIG. 18: The same as the previous Figure, for less central 
events. 



each fluctuating event; and also ii) enhance high-pj^ di- 
rect photon yields due to the high temperature in the 
blobs. 

Since our interest here is just to show the effects of 
fluctuating IC, we chose the simplest decoupling crite- 
rion, namely the usual Cooper-Frye sudden decoupling 
with freezeout temperature Tfo in this Subsection. 
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1. Rapidity and rriT distributions 

Let US first consider Pb + Pb collisions at SPS. In 
Fig. 1161 we show the rapidity distributions for negative 
particles and p — p , respectively, for the most central 
Pb -I- Pb collisions at -^/s = 17.3 A GeV. Each event, com- 
puted from randomly generated IC like the one shown in 
Figs, n and 121 (left), is represented by a thin curve (for 
each type of particle, either negatives or p—p). The thick 
solid lines represent the averages over 50 events, with the 
corresponding dispersions. We compare the average dis- 
tributions with the distributions computed starting from 
the averaged IC like the one shown in Figs. ^ and |21 
(right), represented by dashed lines here. Here, we took 
Tfo = 140 MeV. The data points are shown for compari- 
son [6l |. 

As is seen in this Figure, the rapidity distributions 
show large fluctuations from event to event and, although 
similar, there is a non-negligible difference between the 
average distributions and the ones obtained from the av- 
erage IC, especially in the case of negative particles, con- 
stituted mostly of pions. For the same average initial 
energy, the multiplicity decreases about 6 ~ 7% if fluc- 
tuations exist, confirming what we showed in Sec. 
Finally, the data points in Fig. E| are closer to the aver- 
aged distribution over fluctuating IC. 

What happens with other centralitics is similar as seen 
in Figs. 1171 and 1181 Some differences are: as the colli- 
sions become more peripheral, naturally the dispersions 
increase; ii) the differences between the average distri- 
butions and the ones produced by average IC seems to 
decrease. 

In Fig. 1191 we show the corresponding rriT distri- 
butions, for the most central Pb -I- Pb collisions at 
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FIG. 19: Transverse mass distributions as indicated, com- 
puted at Tfo ~ 140 MeV. The solid lines represent the av- 
erages over the fluctuating distributions (with dispersions), 
whereas the dashed lines are results with the averaged initial 
conditions. The data are from NA49 \6^. 
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FIG. 20: Centrality dependence of transverse mass distribu- FIG. 21: The same as the previous Figure, for less central 

tions, for Pb + Pb collisions at -^s — 17.3 A GeV (SPS). Here, events, 
"with fluctuation" (the solid lines) means results with NeXuS 
fluctuating IC and "without fluctuation" (dashed lines) those 
with the averaged IC. 

with the large fluctuations we saw above of rapidity dis- 
tributions. In part this is due to the logarithmic scale 
y/s — 17.3 A GeV. One can see that apparently the flue- used here. One can conclude, however, that the IC fluc- 
tuation effects on the transverse-momentum spectra are tuations affect very little the slope of the tot spectra, 
very sirrall and the averaged spectra are in good agree- This conclusion is valid for other centralities, as seen in 
ment with those computed with averaged IC, and also Figs. [201 and |^ except for the most peripheral case, 
with data. This smallness of the fluctuation effects on where the dispersions are very large. If we look more 
the transverse-momentum spectra is not in contradiction carefully, we can perceive a small difference between the 
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averaged spectra and those computed with the averaged 
IC, which appears in all the centralities, that is the tail 
of the averaged spectra is more concave and the distri- 
butions become higher, probably because the expansion 
is more violent in this case, due to the high density spots 
in the IC. 

At RHIC energies the results are similar. In Fig. 1221 
we show the pt distributions of charged particles in the 
most central Au + Au collisions at 130 A GeV, calculated 
both with NeXuS fluctuating IC and with averaged IC. 
Like in the case of Pb + Pb collisions at SPS, we used 
sudden freezeout, but with lower freezeout temperature, 
Tfo = 128 MeV, the same value used in (s^l to fit the 
inverse slope parameter of kaons, in the same collisions, 
as shown in Fig. 1151 and Table HI of Section As is seen, 
the fluctuation effects are small in px distribution, both 
curves agree each other and with data, and the aver- 
age over fluctuating events gives a slightly more concave 
shape, just like in Pb -I- Pb collisions at SPS. 

We show in Fig. the pseudorapidity distribution 
of charged particles in the most central Au -I- Au colli- 
sions at 130A GeV, calculated with NeXuS fluctuating IC 
and the corresponding averaged IC, and compared with 
data ^6^. Qualitatively, the behavior is similar to the 
case of Pb -t- Pb collisions at SPS, namely the average f] 
distribution is close to the distribution with the averaged 
IC, being the latter slightly higher than the former. The 
difference here is smaller than at the lower energy case. 
Probably this is due to the decrease of AEi/ < E > in 
Eq. lO as the incident energy increases. 
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FIG. 23: Charged-particle pseudo-rapidity distributions for 
the most central Au-|-Au at 130A GeV. The data are from 
PHOBOS El. 



bution l6£ 



dN 



cos [n{(j) — Tp)] 



(118) 



Thus, 



2. ElUpUc-flow parameter V2 

Let us turn to the elliptic-flow parameter V2 , deflned 
as the second Fourier coefficient of the azimuthal distri- 



Chai-ged Particles; AuAu (130A GeV) 




Pt (GeV) 



FIG. 22: Pt distributions for the most central Au-|-Au at 
130A GeV. The data are from STAR 



W2 =<cos[2(?i-V')] >, (119) 

where the bracket denotes the average value and ^ gives 
the event-dependent collision plane. 

Usually, V2 parameter is interpreted as indicating a 
flow asymmetry caused by the initial-condition asymme- 
try associated with the non-zero impact parameter. In 
this the produced matter in the collision is likely 

to be flattened in the impact-parameter direction, the 
pressure gradient would be larger in this direction and so 
would be the flow. However, as mentioned at the top of 
this section, it was shown in Ref. fl^ that, even for cen- 
tral collisions, fluctuating IC develop azimuthally asym- 
metric flows, because in this case there is no symmetry 
in each event and, experimentally, the impact parameter 
cannot be determined unambiguously. So, one expects 
larger V2 for fluctuating IC as compared to averaged IC 
case. Remark that the fluctuation we are talking about 
is not the often discussed [HE El| flnite- multiplicity fluc- 
tuation at the end of the process. In Fig. |31 results of 
NeXuS-f SPheRIO for V2 of pious produced in Pb-fPb 
collisions at 17. 3A GeV are shown, both with fluctuating 
IC and averaged IC and compared with data [g^ • Large 
discrepancy is seen between the two ways of calculating 
this parameter, and the data are closer to < >, calcu- 
lated with fluctuating IC. 
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Centrality Dependence of (pions; Pb+Pb; 17. 3 A GeV) 
0.08 1 — • — t — • — t — • — t — • — t — • — t — • — t — • — I 



0.07 
0.06 
0.05 
0.04 
0.03 
0.02 
0.01 








o with fluctuation 

■ NA49data 

A without fluctuation 



Centrality 



Preliminary 



correlation function is expressed in terms of the distribu- 
tion function f(x,p) as 



\Ii<l,P)\' 



/(0,pi)/(0,P2) 



(120) 



where P = (pi + P2)/2 and q = (pi ^ P2) and pi is the 
momentum of the zth pion. Usually 

I{q,P) = {alap,)^ f da^P^/(x, P)e^«^ (121) 

In SPH representation, we write /(g, P) as 

^('^' ^) = E 71?^ e-''^?/(-,.P-) , (122) 
bj I iijfj, a J I 



FIG. 24: Centrality dependence ov V2 for pions for Pb+Pb 
collisions at SPS. The data are from NA49 |6l |. 



3. Two-pion interferometry 

As is well known, the identical-particle correlation, also 
known as Hanbury-Brown-Twiss effect (HBT effect) j^^l 
is a powerful tool for probing geometrical sizes of the 
space-time region from which they were emitted. If the 
source is static like a star, it is directly related to the spa- 
tial dimensions of the particle emission source. When ap- 
plied to a dynamical source, however, several non-trivial 
effects appear |6^ |6^, reflecting its time evolution as 
happens in high-energy heavy-ion collisions. Being so, 
the inclusion of IC fluctuations may affect considerably 
the so-called HBT radii, because, as discussed in Sec. El 
the IC in the event-by-event base often show small high- 
density spots in the energy distribution, and our expec- 
tation is that such spots manifest themselves at the end 
when particles are emitted, giving smaller HBT radii. 

We shall discuss here only the recent application of 
NeXuS+SPhcRIO code for HBT effect with IC fluctua- 
tions, for Au+Au collisions at ISOAGeV [gJ. A more 
detailed account of two-particle correlations is given by 
Padula |70l | in a separate paper in this issue. 

For studying the fluctuation effects of IC on HBT cor- 
relation, first we assume Cooper-Frye sudden freezeout 
(FO) at T/.o. = 128 MeV. This freezeout temperature is 
the same one previously found by studying the energy 
dependence of kaon slope parameter T* |5J|, discussed 
in Subsection IV Bl and appears in Table H We showed m 
Subsection IVI A H that indeed this choice of T/.o. repro- 
duces the pT-distribution data at 130A GeV, both with 
averaged and fluctuating IC. We also neglect the reso- 
nance decays. It is argued (tJ that, since resonance 
decays contribute to the correlations with very small q 
values (g <, qmm , where qmin is the minimum measure- 
able g), the experimentally determined HBT radii are 
essentially due to the direct pions. Then the two-particle 



where the sumation is over all the SPH particles. In the 
Cooper-Frye freezeout, these particles should be taken 
where they cross the hyper-surface T = Tf,o. and 
is the normal to this hyper-surface. Notice that, if we 
put Pi ~ P2 — P-, so that g = 0, Eos. 11211 and 11221 are 
reduced, respectively, to Eas. lllOl and lllll that is, to the 
inclusive one-particle distribution. 

In Fig. 1251 we compare the correlation function C2 av- 
eraged over 15 fluctuating events with those computed 
starting from the averaged IC (so, without fluctuations). 
One can see that the IC fluctuations are reflected in large 
fluctuations also in the HBT correlations. When aver- 
aged, the resulting correlation < C2 > are broader than 
those computed with averaged IC, so giving smaller radii 
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FIG. 25: Correlation functions from fluctuating IC and aver- 
aged IC. Sudden freeze-out is used here. The rapidity range 
is —0.5 < Y < 0.5 and qo.a.i which do not appear in the 
horizontal axis are integrated over < qo,s,i < 35MeV. 
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as expected. Also the shape of the correlation functions 
changes. 

We plot the rriT dependence of HBT radii, with Gaus- 
sian fit of C2 , in Fig.Ea together with RHIC data (zllli 
and results with CE, which will be discussed in Sec. 
IVI B 21 It is seen that the smooth IC with sudden FO 
makes the mx dependence of Ro flat or even increasing, 
which is in agreement with other hydro calculations Q 
but in conflict with the data. The fluctuating IC make 
the radii smaller, especially in the case of Ro , however 
without changing the TOy-dependence. 
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FIG. 26; HBT radii and the ratio Ro/ Rs for sudden freeze- 
out (FO) and CE. 1 stands for averaged IC and 2 fluctuating 
IC. Data are from |72. .Ta] . 



to some escaping probability V{x,p) given by Eq. (|114|l . 
There are several nice predictions of the Continuous 
Emission Model (CEM) as discussed in However, 
although more realistic, this description is not handy be- 
cause, V depends on the momentum of the escaping par- 
ticle and, moreover, on the future of the fluid as seen in 
Eq. 1114(1 . In order to make the computation practicable, 
in ,61j] we first took V on the average, i.e., 

r{x,p) ^< r{x,p) >= r{x) . 

Then, approximated linearly the density 
p{x') — as{x') 

(where, for example in the ideal massless pion gas case, 
a = (45C(3))/27r'^ = .278) in the integral of Eq. ^T^. 
Thus, 



(123) 



V{x,p) Vix) — exp ~K 



\ds/dT\ 



where < av > has been included in k = 0.5 a < av >. 
Although approximately, now we can compute V in each 
space-time point. 

In Fig l27l we show the time evolution of the escaping 
probability V{x) given by Ea. H123|l in the mid-rapidity 
plane for the most central Au-|-Au collisions at 130 A 
GeV. The IC have been computed at t = 1 fm and aver- 
aged over 30 NeXuS events. Here and in the computa- 
tions of observables below, the parameter k has been es- 
timated to be .3, corresponding to < av >= 2 fm^ in the 
zero temperature limit and some 20% larger at T = m^r . 
As seen, the probability remains 0.1 < P{x) < 0.8 in a 
quite large domain, especially for t > 10 fm, indicating 
that both the emission zone and duration are expected 
to be large, in opposition to the standard sudden freeze- 
out case. For comparison, we show also the temperature 
distribution in Fig l27l with some isotherms. 

To calculating the spectra, now Eq. H117|l is trans- 
lated into SPH language, which is given by Eq. (|111|1 . 
but where the sum is computed in the present case not 
over T = T/.o. hypersurface but picking out SPH particles 
according to the probability V, given by Eq. H123|) with 
the normal pointing to the 4-gradient of V. Since 
our procedure favors emission from fast outgoing SPH 
"particles" , because p decreases faster there and so does 
s in this case making V larger, we believe that the main 
feature of CEM is preserved in our approximation. 



niT distributions 



B. Effects of continuous emission 

As discussed in Sec. it is more likely that the de- 
coupling occurs not suddenly, but continuously from any 
point of the fluid and at any instant of time, according 



In Fig. I2HI we show the charged mr distribution at 
mid-rapidity in the most central Au + Au collisions at 
130A GeV, computed by using CEM, with and without 
fluctuations. Compare this figure with Fig. 1221 where the 
same data are compared with calculations using sudden 
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i;(fm) 

FIG. 27: Upper panel: Probability distribution as given by 
eg. 111231 for the most central Au+Au collisions at 130 A GeV, 
in the mid-rapidity plane for averaged IC. Lower panel: Cor- 
responding temperature distribution. 



freezeout. One sees that, although the emission mecha- 
nisms are quite different, the results are similar, for the 
choice of the parameters, T/o in one case and k in the 
other. However, the origins of the spectrum shape are 
different. Whereas, in the freezeout case, all the parti- 
cles are emitted at the same temperature Tfo and the 




Pt (GeV) 

FIG. 28: Charged-particle pr distributions in GEM for the 
most central Au + Au at 130A GeV. The data are from 
STAR |6^. 



concave shape of the spectra is due to different trans- 
verse velocities of the fluid at different instants of time, 
in the continuous emission case, the large-pr particles 
are in general emitted earlier and at higher temperature, 
the small-pT particles are emitted later, when the fluid 
is cooler. 



2. Two-pion inter ferometry 

Let us now consider the effect of continuous emission 
on HBT effect. Since HBT is sensitive to the space- 
time geometry of the fluid, and GEM produces impor- 
tant modiflcation of the emission zone, we expect consid- 
erable changes in the so-called HBT radii. As mentioned 
above, according to this picture, the large-pT particles 
are mainly emitted at early times when the fluid is hot 
and mostly from its surface, whereas the small-pr com- 
ponents are emitted later when the fluid is cooler and 
from larger spatial domain. 

In |23|, we considered this effect, using the boost- 
invariant solution [30| . and showed that whereas the so- 
called side radius is independent of the average pt , the 
out radius decreases with < px >, because of the reason 
mentioned above. This behavior is expected to essentially 
remain in the general 3-dimensional expansion, described 
by SPheRIO code. 

To compute the correlation function C2(q, P) in GEM, 
we flrst rewrite the integral H121|l as 

J (70 

+ J d^xd^[P^ffree{x,PW'^^, (124) 

which is similar to Eq. (|117|) . Then, translate it into SPH 
language by using Eq. H122() , but with the sum evaluated 
by picking out SPH particles according to the probability 
V, given by Eq. H123(l with the normal n^^ pointing to 
the 4-gradient of V, exactly in the same way as done in 
calculating the inclusive spectrum. 

We plot in Fig. 1201 some results for the dependence 
of HBT radii computed in this way, both with fluctuating 
IC and averaged IC, for Au -I- Au collisions at 130A GeV. 
For comparison, we show also the results with sudden 
freezeout and data [llllll. Comparing the averaged IG 
case (with GEM (GEl)), with the corresponding freeze- 
out (FOl), one sees that, while Rl remains essentially 
the same, Rg decreases faster and as for Rg , it decreases 
now inverting its behavior. 

In Fig.|2ni one can also see the combined effect of fluc- 
tuating initial conditions together with continuous emis- 
sion (curve GE2). All the radii are smaller than CEl, 
with averaged IG (but using GEM), as happend in the 
sudden freezeout case. The agreement with data is exce- 
lent both for i?^ and Rs , and improved considerably the 
results for Ro with respect to the usual hydrodynamic 
description. 
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VII. CONCLUSIONS AND OUTLOOK 

In this survey, we discussed on several aspects of ap- 
plications of hydrodynamic model to nucleus-nucleus col- 
lisons, giving especial emphases on i) method of solving 
the hydrodynamic equations for arbitrary configurations; 
ii) accounting for the probable event-by-event fluctua- 
tions of the initial conditions; and iii) the decoupling cri- 
teria for obtaining the observables. 

The Smoothed-Particle Hydrodynamic approach, us- 
ing a special hyperbolic coordinate system, is particularly 
interesting tool for describing rapidly expandig matter 
as that in high-energy nucleus-nucleus coUisons, because 
of its efficiency and flexibility, as demonstrated through 
examples considered in Sec. IIV C 51 IIV C 61 and through 
applications, shown in Sec. IVII where the systems, in 
general, do not present any symmetry. 

The initial conditions for hydrodynamic expansion 
of matter formed in relativistic heavy-ion collisions are 
likely neither smooth nor symmetrical, because of small 
size of these systems. This property has also been shown 
by some event generators, which take the microscopic 
dynamics into account. The fluctuations in the initial 
conditions may produce large discrepancies in the com- 
puted results for some observables, in comparison with 
those obtained with the usual averaged, smooth and sym- 
metrical initial conditions. So, it is our opinion that, to 
extracting correct information from certain experimental 
data, the inclusion of the fluctuations mentioned above 
is mandatory. 

We find that, although Cooper- Frye formula can give 
good results in many cases, a more realistic decoupling 
description is needed for obtaining many other observ- 
ables. A description we proposed is the Continuous Emis- 
sion Model. One typical example we showed here, in 
which this model improve considerably the description 
of data is the my dependence of HBT radii, observed at 
RHIC. 

In addition to the points summarized above, there are 
several basic questions to be addressed in the hydrody- 
namical scenary. First of all, although for the bulk prop- 
erties of observables so far studied the ideal hydrodynam- 
ical description works fairly well, we still do not know 
how the non-equilibrium processes affect the results of 
these studies. This question refers both to the initial 
conditions and to the final particle decoupling process. 
For example, the present study of continuous emission 
showed that the real particle emission mechanism seems 
to be very far from the conventional Cooper- Frye scenary 
of sharp freeze-out surface. In Fig.|23 we showed that the 
escape probability V varies rather slowly both in space 
and time, and there is no indication of a sharp freezeout 
hypersurface. Another way of showing the same is given 
in Fig. 1291 where we plotted the entropy vs. temperature 
of the SPH particles at the emission point. The tem- 
perature values corresponding to the particle emission 
spread widely. This shows that, even for a large system 
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FIG. 29: The entropy liberated vs. temperature of SPH par- 
ticles at the decoupling point. The freeze-out temperature 
T/o used in the studies above is also shown. 



(Au-|-Au) there exists a substantial part of the system 
which remains out of equilibrium for a long time. Up 
to now, no dynamical reaction of these non-equilibrium 
components on the hydrodynamical evolution is studied. 
In this sense, it is very irnpqrtant to develop transport 
theoretical investigations [t^ [t^ related to hydro de- 
scription. 

Another interesting type of phenomena related to the 
hydrodynamical model is the ocurrence of instabilities 
associated to the phase transition and surface tension 
(finite-size effects), such as spinodal instabilities |77j. 
Some studies in this direction is in progress. 
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